"""
美元期货指数聚类分析脚本

本脚本依据《任务要求.md》实现三种聚类方案：
1) 直接对标准化特征进行 KMeans 聚类；
2) PCA 降维后进行 KMeans 聚类；
3) 因子分析（最大似然 + Varimax旋转）降维后进行 KMeans 聚类。

功能特性：
- 数据清洗与预处理：日期解析、交易量单位转换、涨跌幅百分号处理；
- 特征选择与标准化：收盘/开盘/高/低/交易量/涨跌幅；
- K 值选择：在 K=2..10 范围内以轮廓系数选择最优；
- 评估与对比：轮廓系数（样本级与簇级）、时间复杂度（运行时间）、稳定性（多次运行方差）、业务可解释性评分；
- 可视化：轮廓系数折线图、PCA累计方差（Scree）图、主成分载荷热力图、因子载荷热力图、二维散点图（PCA/FA 得分空间）；
- 报告生成：Markdown 格式，包含公式、代码片段、对比表与结论推荐；
- Mac中文显示支持：自动设置中文字体（PingFang SC/SimHei）与负号显示。

用法示例：
python3 homework/chapter6/cluster_analysis.py \
  --csv homework/chapter6/美元指数期货历史数据.csv \
  --out-dir homework/chapter6/outputs \
  --k-min 2 --k-max 10 --random-state 123 --n-init 10 --var-threshold 0.7

"""

from __future__ import annotations

import argparse
import os
import time
import math
from dataclasses import dataclass
from typing import Dict, List, Tuple, Optional

import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, silhouette_samples
from sklearn.decomposition import PCA, FactorAnalysis


# =====================
# 工具与可视化设置
# =====================

def set_chinese_font() -> None:
    """
    设置 Matplotlib 中文字体，优先使用 Mac 上常见的 'PingFang SC'，
    若不可用则回退为 'SimHei'。并修正负号显示。
    """
    try:
        matplotlib.rcParams["axes.unicode_minus"] = False
        # 优先使用苹方（Mac内置），否则回退黑体
        font_candidates = ["PingFang SC", "Songti SC", "Heiti SC", "SimHei"]
        from matplotlib.font_manager import fontManager

        available = {f.name for f in fontManager.ttflist}
        for name in font_candidates:
            if name in available:
                matplotlib.rcParams["font.sans-serif"] = [name]
                return
        # 如果都不可用，仍然设置一个常见中文字体名以尽可能显示中文
        matplotlib.rcParams["font.sans-serif"] = ["SimHei"]
    except Exception:
        # 安静失败，不影响后续绘图
        matplotlib.rcParams["font.sans-serif"] = ["SimHei"]


def ensure_dirs(out_dir: str) -> Dict[str, str]:
    """
    创建输出目录结构，包括 figures 子目录。

    参数：
    - out_dir: 输出根目录。

    返回：
    - 路径字典，包含 'out_dir' 与 'figures_dir'。
    """
    figures_dir = os.path.join(out_dir, "figures")
    os.makedirs(figures_dir, exist_ok=True)
    return {"out_dir": out_dir, "figures_dir": figures_dir}


# =====================
# 数据处理
# =====================

def _parse_chinese_date(s: str) -> pd.Timestamp:
    """
    将形如 '2024年5月31日' 的中文日期解析为 Pandas 时间戳。
    如果解析失败，回退为通用 pd.to_datetime。
    """
    try:
        # 简单替换，避免locale依赖
        s = s.replace("年", "-").replace("月", "-").replace("日", "")
        return pd.to_datetime(s, format="%Y-%m-%d")
    except Exception:
        return pd.to_datetime(s, errors="coerce")


def _parse_volume(s: str) -> float:
    """
    交易量字符串转为数值。
    支持如 '13.90K'（千）与 '1.2M'（百万），以及纯数字。
    """
    if s is None or (isinstance(s, float) and math.isnan(s)):
        return np.nan
    s = str(s).strip()
    try:
        if s.endswith("K") or s.endswith("k"):
            return float(s[:-1]) * 1_000.0
        if s.endswith("M") or s.endswith("m"):
            return float(s[:-1]) * 1_000_000.0
        return float(s)
    except Exception:
        return np.nan


def _parse_pct(s: str) -> float:
    """
    涨跌幅百分数字符串转小数。例如 '-0.03%' -> -0.0003。
    """
    if s is None or (isinstance(s, float) and math.isnan(s)):
        return np.nan
    s = str(s).strip().replace("%", "")
    try:
        return float(s) / 100.0
    except Exception:
        return np.nan


def load_and_preprocess(csv_path: str) -> pd.DataFrame:
    """
    读取并清洗数据：解析日期、交易量与涨跌幅，按日期升序排序。

    参数：
    - csv_path: CSV文件路径。

    返回：
    - 清洗后的 DataFrame，包含原始列与标准化前的数值列。
    """
    df = pd.read_csv(csv_path, encoding="utf-8")
    # 列名期望：日期、收盘、开盘、高、低、交易量、涨跌幅
    df["日期"] = df["日期"].apply(_parse_chinese_date)
    df["交易量"] = df["交易量"].apply(_parse_volume)
    df["涨跌幅"] = df["涨跌幅"].apply(_parse_pct)
    # 转浮点
    for col in ["收盘", "开盘", "高", "低"]:
        df[col] = pd.to_numeric(df[col], errors="coerce")
    # 去除缺失并排序
    df = df.dropna(subset=["日期"]).sort_values("日期").reset_index(drop=True)
    return df


def select_numeric_features(df: pd.DataFrame) -> pd.DataFrame:
    """
    选择用于建模的数值型特征列。
    包括：收盘、开盘、最高、最低、交易量、涨跌幅。
    """
    cols = ["收盘", "开盘", "高", "低", "交易量", "涨跌幅"]
    return df[cols].copy()


def standardize_features(X: pd.DataFrame) -> Tuple[np.ndarray, StandardScaler]:
    """
    使用 StandardScaler 对特征进行标准化。

    返回：
    - X_scaled: 标准化后的特征矩阵 (numpy.ndarray)
    - scaler: 拟合的 StandardScaler，用于逆变换解释簇中心。
    """
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X.values)
    return X_scaled, scaler


# =====================
# 评估与可视化
# =====================

def silhouette_k_search(
    X: np.ndarray,
    k_min: int,
    k_max: int,
    random_state: int,
    n_init: int,
) -> Tuple[int, List[int], List[float]]:
    """
    在给定 K 范围内搜索最佳簇数（按平均轮廓系数最大）。

    返回：
    - best_k: 得分最高的 K 值
    - ks: 测试的 K 列表
    - scores: 对应的轮廓系数列表
    """
    ks, scores = [], []
    for k in range(max(2, k_min), k_max + 1):
        km = KMeans(n_clusters=k, random_state=random_state, n_init=n_init)
        labels = km.fit_predict(X)
        score = silhouette_score(X, labels)
        ks.append(k)
        scores.append(float(score))
    best_idx = int(np.argmax(scores))
    return ks[best_idx], ks, scores


def plot_silhouette_curve(
    ks: List[int],
    scores: List[float],
    title: str,
    out_path: str,
) -> None:
    """
    绘制轮廓系数随 K 值变化曲线。
    输出 PNG（dpi≥300）。
    """
    set_chinese_font()
    plt.figure(figsize=(7, 4), dpi=300)
    plt.plot(ks, scores, marker="o")
    plt.title(title)
    plt.xlabel("簇数 K")
    plt.ylabel("平均轮廓系数")
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig(out_path)
    plt.close()


def silhouette_details(X: np.ndarray, labels: np.ndarray) -> Tuple[float, Dict[int, float], np.ndarray]:
    """
    计算轮廓系数的整体、簇级平均与样本级细节。
    """
    overall = silhouette_score(X, labels)
    sample_scores = silhouette_samples(X, labels)
    cluster_scores: Dict[int, float] = {}
    for c in np.unique(labels):
        cluster_scores[int(c)] = float(np.mean(sample_scores[labels == c]))
    return float(overall), cluster_scores, sample_scores


def pca_reduce(
    X: np.ndarray,
    feature_names: List[str],
    var_threshold: float = 0.7,
) -> Tuple[np.ndarray, PCA, np.ndarray, np.ndarray]:
    """
    PCA 降维，保留累计方差贡献率≥阈值的主成分。

    返回：
    - Z: 主成分得分矩阵
    - pca: 拟合的 PCA 对象
    - explained_ratio: 每个主成分的方差贡献率
    - loadings: 主成分载荷矩阵（形状：n_features × n_components）
    """
    # 先拟合完整 PCA 以确定主成分数
    full = PCA()
    full.fit(X)
    cum = np.cumsum(full.explained_variance_ratio_)
    n_components = int(np.searchsorted(cum, var_threshold) + 1)
    pca = PCA(n_components=n_components)
    Z = pca.fit_transform(X)
    explained_ratio = pca.explained_variance_ratio_
    # 载荷：components_ 形状为 (n_components × n_features)，转置为 (n_features × n_components)
    loadings = pca.components_.T
    return Z, pca, explained_ratio, loadings


def plot_pca_scree(
    explained_ratio_full: np.ndarray,
    title: str,
    out_path: str,
) -> None:
    """
    绘制 PCA Scree 图（单个主成分贡献与累计贡献）。
    """
    set_chinese_font()
    plt.figure(figsize=(7, 4), dpi=300)
    x = np.arange(1, len(explained_ratio_full) + 1)
    plt.bar(x, explained_ratio_full, label="单个主成分贡献率")
    plt.plot(x, np.cumsum(explained_ratio_full), marker="o", label="累计贡献率")
    plt.xlabel("主成分编号")
    plt.ylabel("方差贡献率")
    plt.title(title)
    plt.legend()
    plt.tight_layout()
    plt.savefig(out_path)
    plt.close()


def plot_heatmap(
    matrix: np.ndarray,
    row_labels: List[str],
    col_labels: List[str],
    title: str,
    out_path: str,
) -> None:
    """
    绘制载荷矩阵热力图。
    """
    set_chinese_font()
    plt.figure(figsize=(6, 4), dpi=300)
    sns.heatmap(
        pd.DataFrame(matrix, index=row_labels, columns=col_labels),
        cmap="RdBu_r",
        center=0.0,
        annot=True,
        fmt=".2f",
    )
    plt.title(title)
    plt.tight_layout()
    plt.savefig(out_path)
    plt.close()


def plot_scatter_2d(
    Z: np.ndarray,
    labels: np.ndarray,
    title: str,
    out_path: str,
    axis_labels: Tuple[str, str] = ("维度1", "维度2"),
) -> None:
    """
    在二维得分空间绘制散点图，并按簇上色。
    """
    set_chinese_font()
    plt.figure(figsize=(6, 5), dpi=300)
    palette = sns.color_palette("tab10", n_colors=len(np.unique(labels)))
    for c in np.unique(labels):
        mask = labels == c
        plt.scatter(Z[mask, 0], Z[mask, 1], s=25, label=f"簇 {int(c)}", color=palette[int(c) % len(palette)])
    plt.title(title)
    plt.xlabel(axis_labels[0])
    plt.ylabel(axis_labels[1])
    plt.legend()
    plt.tight_layout()
    plt.savefig(out_path)
    plt.close()


# =====================
# 因子分析与旋转
# =====================

def varimax(Phi: np.ndarray, gamma: float = 1.0, q: int = 20, tol: float = 1e-6) -> np.ndarray:
    """
    Varimax 旋转实现（正交旋转），用于提高因子可解释性。

    参数：
    - Phi: 载荷矩阵（形状：n_features × n_factors）
    - gamma: 旋转参数（默认1.0）
    - q: 最大迭代次数
    - tol: 收敛阈值

    返回：
    - 旋转后的载荷矩阵（与 Phi 同形状）
    """
    p, k = Phi.shape
    R = np.eye(k)
    d = 0
    for _ in range(q):
        d_old = d
        Lambda = np.dot(Phi, R)
        u, s, vh = np.linalg.svd(
            np.dot(Phi.T, Lambda**3 - (gamma / p) * np.dot(Lambda, np.diag(np.diag(np.dot(Lambda.T, Lambda)))))
        )
        R = np.dot(u, vh)
        d = np.sum(s)
        if d_old != 0 and d / d_old < 1 + tol:
            break
    return np.dot(Phi, R)


def choose_n_factors_by_corr_eigen(X: np.ndarray, threshold: float = 0.7) -> int:
    """
    依据相关矩阵的特征值累计贡献率选择因子数，使累计贡献度≥阈值。
    """
    C = np.corrcoef(X, rowvar=False)
    evals, _ = np.linalg.eigh(C)
    evals = np.flip(np.sort(evals))  # 从大到小
    cum = np.cumsum(evals) / np.sum(evals)
    n = int(np.searchsorted(cum, threshold) + 1)
    return max(1, n)


def factor_analysis_reduce(
    X: np.ndarray,
    feature_names: List[str],
    var_threshold: float,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, int, float]:
    """
    因子分析降维流程：
    1. 基于相关矩阵特征值选择因子数量（累计解释度≥阈值）
    2. 最大似然估计
    3. Varimax 旋转载荷矩阵
    
    返回：
    - Z: 因子得分（n_samples × n_factors）
    - loadings_rot: 旋转后因子载荷
    - loadings_raw: 原始因子载荷（未旋转）
    - n_factors: 因子数量
    - approx_explained: 近似累计解释度（来自相关矩阵特征值法）
    """
    n_factors = choose_n_factors_by_corr_eigen(X, threshold=var_threshold)
    fa = FactorAnalysis(n_components=n_factors)  # 移除rotation参数
    Z = fa.fit_transform(X)
    # sklearn 的因子载荷可通过 components_ 获得（n_components × n_features），转置为 n_features × n_components
    loadings_raw = fa.components_.T
    loadings_rot = varimax(loadings_raw)
    # 近似累计解释度（用于报告）：
    C = np.corrcoef(X, rowvar=False)
    evals, _ = np.linalg.eigh(C)
    evals = np.flip(np.sort(evals))
    approx_explained = float(np.sum(evals[:n_factors]) / np.sum(evals))
    return Z, loadings_rot, loadings_raw, n_factors, approx_explained


# =====================
# 业务解释与评分
# =====================

def label_clusters_by_centers(
    centers_z: np.ndarray,
    feature_names: List[str],
) -> Dict[int, str]:
    """
    根据簇中心的标准化值为簇命名，提供业务解释标签。
    简化规则：
    - 涨跌幅高正 → “上涨期”；低负 → “下跌期”；绝对值大 → “高波动期”；
    - 交易量高 → “高交易量”；低 → “低交易量”；
    - 收盘价高 → “高价位”；低 → “低价位”。
    """
    idx_price = feature_names.index("收盘")
    idx_change = feature_names.index("涨跌幅")
    idx_volume = feature_names.index("交易量")
    labels: Dict[int, str] = {}
    for i, z in enumerate(centers_z):
        parts: List[str] = []
        # 涨跌方向/波动
        if z[idx_change] >= 0.3:
            parts.append("上涨期")
        elif z[idx_change] <= -0.3:
            parts.append("下跌期")
        else:
            parts.append("震荡期")
        # 波动强弱（以涨跌幅绝对值衡量）
        if abs(z[idx_change]) >= 0.6:
            parts.append("高波动")
        elif abs(z[idx_change]) <= 0.2:
            parts.append("低波动")
        # 交易量
        if z[idx_volume] >= 0.5:
            parts.append("高交易量")
        elif z[idx_volume] <= -0.5:
            parts.append("低交易量")
        # 价位
        if z[idx_price] >= 0.5:
            parts.append("高价位")
        elif z[idx_price] <= -0.5:
            parts.append("低价位")
        labels[i] = "、".join(parts) if parts else f"簇 {i}"
    return labels


def top_loadings(matrix: np.ndarray, feature_names: List[str], top_n: int = 3) -> List[List[Tuple[str, float]]]:
    """
    提取每个维度（主成分或因子）的载荷TopN（按绝对值排序）。
    返回：列表，长度为维度数，每项为 [(特征名, 载荷值), ...]。
    """
    result: List[List[Tuple[str, float]]] = []
    for j in range(matrix.shape[1]):
        col = matrix[:, j]
        idx = np.argsort(np.abs(col))[::-1][:top_n]
        result.append([(feature_names[i], float(col[i])) for i in idx])
    return result


def name_dimension_by_tops(tops: List[Tuple[str, float]], is_factor: bool = False) -> Tuple[str, str]:
    """
    根据Top载荷特征给主成分/因子命名，并返回简短业务解释。
    规则：
    - 含“交易量”居前 → 成交量主成分/因子（市场活跃度）
    - 含“涨跌幅”居前 → 动量/波动主成分/因子（价格变动强度）
    - 多个价格相关（收盘/开盘/高/低）居前且同号 → 价格水平主成分/因子（整体价格水准）
    - 价格相关为混合符号 → 价格振幅/形态主成分/因子（高低价差异）
    """
    names_prices = {"收盘", "开盘", "高", "低"}
    top_feats = [t[0] for t in tops]
    top_vals = [t[1] for t in tops]
    is_comp = not is_factor
    prefix = "主成分" if is_comp else "因子"

    if "交易量" in top_feats:
        return f"成交量{prefix}", "以交易量为主导，刻画市场活跃度与参与度"
    if "涨跌幅" in top_feats:
        return f"动量/波动{prefix}", "以涨跌幅为主导，反映价格变动方向与强度"

    price_mask = [1 if f in names_prices else 0 for f in top_feats]
    if sum(price_mask) >= 2:
        # 判断符号一致性
        same_sign = np.all(np.sign(top_vals) == np.sign(top_vals[0]))
        if same_sign:
            return f"价格水平{prefix}", "多个价格维度同向载荷，反映整体价格水准"
        else:
            return f"价格形态/振幅{prefix}", "价格维度载荷符号不一，反映高低价差与形态"

    # 兜底：按第一特征命名
    return f"{top_feats[0]}{prefix}", "以该特征为主导的综合维度"


def name_all_dimensions(loadings: np.ndarray, feature_names: List[str], is_factor: bool = False) -> List[Tuple[str, List[Tuple[str, float]], str]]:
    """
    为全部主成分/因子生成命名与解释。
    返回：[(名称, Top载荷[(特征,值)], 解释), ...]
    """
    tops = top_loadings(loadings, feature_names, top_n=3)
    named: List[Tuple[str, List[Tuple[str, float]], str]] = []
    for j in range(loadings.shape[1]):
        nm, desc = name_dimension_by_tops(tops[j], is_factor=is_factor)
        named.append((nm, tops[j], desc))
    return named


def interpretability_score_from_loadings(loadings: np.ndarray) -> int:
    """
    基于载荷矩阵的稀疏性与显著性，粗略评估可解释性（1-5分）。
    规则：统计 |loading|>=0.5 的条目数，越多得分越高；过于分散略降权。
    """
    strong = np.sum(np.abs(loadings) >= 0.5)
    if strong >= 8:
        return 5
    if strong >= 5:
        return 4
    if strong >= 3:
        return 3
    if strong >= 1:
        return 2
    return 1


def interpretability_score_kmeans(centers_z: np.ndarray) -> int:
    """
    对直接KMeans的可解释性打分：
    依据簇中心在关键维度（交易量/涨跌幅/收盘）的分化程度。
    """
    key_dims = [0, 4, 5]  # 收盘(0)、交易量(4)、涨跌幅(5) —— 与 feature_names 对应
    spread = np.mean(np.std(centers_z[:, key_dims], axis=0))
    if spread >= 0.8:
        return 5
    if spread >= 0.6:
        return 4
    if spread >= 0.4:
        return 3
    if spread >= 0.2:
        return 2
    return 1


# =====================
# 稳定性评估
# =====================

def stability_variance(
    X: np.ndarray,
    k: int,
    seeds: List[int],
    n_init: int,
) -> float:
    """
    在固定簇数 k 下，以不同随机种子重复运行，计算平均轮廓系数的方差。
    注意：主结果仍以 random_state=123 输出，本函数仅用于稳定性对比。
    """
    scores: List[float] = []
    for rs in seeds:
        km = KMeans(n_clusters=k, random_state=rs, n_init=n_init)
        labels = km.fit_predict(X)
        scores.append(float(silhouette_score(X, labels)))
    return float(np.var(scores, ddof=1)) if len(scores) > 1 else 0.0


# =====================
# 报告生成
# =====================

def write_markdown_report(
    out_dir: str,
    figures: Dict[str, str],
    summary: Dict[str, Dict[str, float]],
    pca_info: Dict[str, any],
    fa_info: Dict[str, any],
    code_params: Dict[str, str],
    cluster_labels_direct: Dict[int, str],
    cluster_labels_pca: Dict[int, str],
    cluster_labels_fa: Dict[int, str],
    pca_named: List[Tuple[str, List[Tuple[str, float]], str]],
    fa_named: List[Tuple[str, List[Tuple[str, float]], str]],
) -> str:
    """
    生成 Markdown 报告文件，包含方法原理、代码片段、图表、对比与结论。
    返回报告文件路径。
    """
    md_path = os.path.join(out_dir, "cluster_analysis_report.md")

    def _fmt(x: float) -> str:
        return f"{x:.4f}"

    pca_explained = pca_info.get("explained_ratio_full", None)
    pca_components = pca_info.get("n_components", None)
    fa_components = fa_info.get("n_factors", None)
    fa_explained = fa_info.get("approx_explained", None)

    lines: List[str] = []
    lines.append("# 美元期货指数聚类分析报告")
    lines.append("")
    lines.append("## 1. 方法原理")
    lines.append("- KMeans：最小化簇内平方误差，核心目标函数为 $\\min \\sum_{i} \\lVert x_i - \\mu_{c(i)} \\rVert^2$。")
    lines.append("- 轮廓系数：$s(i) = \\frac{b(i) - a(i)}{\\max(a(i), b(i))}$，其中 $a(i)$ 为样本到同簇内的平均距离，$b(i)$ 为到最近其他簇的平均距离。")
    lines.append("- PCA：在标准化空间中，令协方差矩阵特征分解得到主成分，选择累计方差贡献率≥阈值的前若干主成分。")
    lines.append("- 因子分析（ML+Varimax）：以最大似然估计因子载荷并进行 Varimax 旋转，使载荷呈块状集中以提升可解释性。")
    lines.append("")
    lines.append("## 2. 关键参数与代码片段")
    lines.append("```python")
    lines.append(f"RANDOM_STATE = {code_params['random_state']}  # 所有主流程固定种子")
    lines.append(f"N_INIT = {code_params['n_init']}            # KMeans多次初始化")
    lines.append("# 直接聚类")
    lines.append("KMeans(n_clusters=K, random_state=RANDOM_STATE, n_init=N_INIT)")
    lines.append("# PCA降维")
    lines.append("PCA(n_components=auto_by_threshold)  # 累计方差≥70%")
    lines.append("# FA降维（最大似然 + Varimax）")
    lines.append("FactorAnalysis(n_components=auto_by_corr_eigen); loadings = varimax(loadings)")
    lines.append("```")
    lines.append("")
    lines.append("## 3. 可视化结果")
    lines.append(f"- 直接聚类轮廓系数曲线：![]({os.path.relpath(figures['kmeans_curve'], out_dir)})")
    lines.append(f"- PCA Scree图（累计方差）：![]({os.path.relpath(figures['pca_scree'], out_dir)})")
    lines.append(f"- PCA载荷热力图：![]({os.path.relpath(figures['pca_loadings'], out_dir)})")
    lines.append(f"- PCA二维散点：![]({os.path.relpath(figures['pca_scatter'], out_dir)})")
    lines.append(f"- PCA聚类轮廓曲线：![]({os.path.relpath(figures['pca_curve'], out_dir)})")
    lines.append(f"- 因子载荷热力图：![]({os.path.relpath(figures['fa_loadings'], out_dir)})")
    lines.append(f"- 因子二维散点：![]({os.path.relpath(figures['fa_scatter'], out_dir)})")
    lines.append(f"- 因子聚类轮廓曲线：![]({os.path.relpath(figures['fa_curve'], out_dir)})")
    lines.append("")
    lines.append("## 4. 评估与业务解释")
    lines.append("- 直接聚类：整体轮廓系数=" + _fmt(summary['kmeans']['silhouette']) + 
                 ", 运行时间(s)=" + _fmt(summary['kmeans']['runtime']) + 
                 ", 稳定性方差=" + _fmt(summary['kmeans']['stability_var']) + 
                 ", 可解释性评分=" + str(int(summary['kmeans']['interpretability'])))
    lines.append("- PCA + KMeans：整体轮廓系数=" + _fmt(summary['pca']['silhouette']) + 
                 ", 运行时间(s)=" + _fmt(summary['pca']['runtime']) + 
                 ", 稳定性方差=" + _fmt(summary['pca']['stability_var']) + 
                 ", 可解释性评分=" + str(int(summary['pca']['interpretability'])) + 
                 f"，主成分数={pca_components}")
    lines.append("- FA + KMeans：整体轮廓系数=" + _fmt(summary['fa']['silhouette']) + 
                 ", 运行时间(s)=" + _fmt(summary['fa']['runtime']) + 
                 ", 稳定性方差=" + _fmt(summary['fa']['stability_var']) + 
                 ", 可解释性评分=" + str(int(summary['fa']['interpretability'])) + 
                 f"，因子数={fa_components}，近似累计解释度={_fmt(fa_explained)}")
    lines.append("")
    lines.append("### 4.1 簇命名与业务解释")
    lines.append("- 直接KMeans：")
    for cid, name in cluster_labels_direct.items():
        lines.append(f"  - 簇 {cid}: {name}")
    lines.append("- PCA+KMeans：")
    for cid, name in cluster_labels_pca.items():
        lines.append(f"  - 簇 {cid}: {name}")
    lines.append("- FA+KMeans：")
    for cid, name in cluster_labels_fa.items():
        lines.append(f"  - 簇 {cid}: {name}")
    lines.append("")
    lines.append("### 4.2 主成分与因子命名")
    lines.append("- PCA 主成分命名：")
    for i, (nm, tops, desc) in enumerate(pca_named, start=1):
        pairs = ", ".join([f"{f}({v:.2f})" for f, v in tops])
        lines.append(f"  - PC{i}: {nm}；主要载荷：{pairs}；解释：{desc}")
    lines.append("- 因子分析 因子命名：")
    for i, (nm, tops, desc) in enumerate(fa_named, start=1):
        pairs = ", ".join([f"{f}({v:.2f})" for f, v in tops])
        lines.append(f"  - F{i}: {nm}；主要载荷：{pairs}；解释：{desc}")
    lines.append("")
    lines.append("## 5. 对比表")
    lines.append("| 方法 | 轮廓系数 | 运行时间(s) | 稳定性方差 | 可解释性评分 |")
    lines.append("|------|---------:|------------:|-----------:|------------:|")
    lines.append(f"| 直接KMeans | {_fmt(summary['kmeans']['silhouette'])} | {_fmt(summary['kmeans']['runtime'])} | {_fmt(summary['kmeans']['stability_var'])} | {int(summary['kmeans']['interpretability'])} |")
    lines.append(f"| PCA+KMeans | {_fmt(summary['pca']['silhouette'])} | {_fmt(summary['pca']['runtime'])} | {_fmt(summary['pca']['stability_var'])} | {int(summary['pca']['interpretability'])} |")
    lines.append(f"| FA+KMeans  | {_fmt(summary['fa']['silhouette'])} | {_fmt(summary['fa']['runtime'])} | {_fmt(summary['fa']['stability_var'])} | {int(summary['fa']['interpretability'])} |")
    lines.append("")
    lines.append("## 6. 结论与推荐")
    lines.append("- 若关注可解释性与维度压缩，建议采用 FA+KMeans（Varimax后载荷更具块状结构）；")
    lines.append("- 若更注重整体聚类分辨率与运行效率，PCA+KMeans 通常表现稳健；")
    lines.append("- 直接KMeans在原空间可作为基线方案，用于与降维后的聚类效果做对照。")
    lines.append("")
    lines.append("> 注：主流程均以 `random_state=123, n_init=10` 固定，稳定性指标为固定 K 下在不同种子重复运行所得的轮廓系数方差，仅用于对比，不影响主结果复现性。")

    with open(md_path, "w", encoding="utf-8") as f:
        f.write("\n".join(lines))
    return md_path


# =====================
# 主流程
# =====================

def run_pipeline(
    csv_path: str,
    out_dir: str,
    k_min: int,
    k_max: int,
    random_state: int,
    n_init: int,
    var_threshold: float,
) -> Tuple[str, Dict[str, str]]:
    """
    执行完整分析流程，生成图表与报告。

    返回：
    - 报告文件路径
    - 生成图表路径字典
    """
    paths = ensure_dirs(out_dir)
    figures_dir = paths["figures_dir"]

    # 1) 数据加载与标准化
    t0 = time.perf_counter()
    df = load_and_preprocess(csv_path)
    X_df = select_numeric_features(df)
    feature_names = list(X_df.columns)
    X_scaled, scaler = standardize_features(X_df)
    # 2) 直接KMeans
    t1 = time.perf_counter()
    k_best_1, ks1, sc1 = silhouette_k_search(X_scaled, k_min, k_max, random_state, n_init)
    curve1 = os.path.join(figures_dir, "kmeans_silhouette_curve.png")
    plot_silhouette_curve(ks1, sc1, "方式一：直接聚类的轮廓系数-簇数关系", curve1)
    km1 = KMeans(n_clusters=k_best_1, random_state=random_state, n_init=n_init)
    labels1 = km1.fit_predict(X_scaled)
    sil1_overall, sil1_cluster, sil1_samples = silhouette_details(X_scaled, labels1)
    centers1_z = km1.cluster_centers_
    clusters1_label = label_clusters_by_centers(centers1_z, feature_names)
    stability1 = stability_variance(X_scaled, k_best_1, [101, 123, 231, 345, 567], n_init)
    t2 = time.perf_counter()

    # 3) PCA 降维后聚类
    t3 = time.perf_counter()
    # Scree基于完整拟合
    pca_full = PCA()
    pca_full.fit(X_scaled)
    scree_path = os.path.join(figures_dir, "pca_scree.png")
    plot_pca_scree(pca_full.explained_variance_ratio_, "PCA累计方差贡献率（Scree图）", scree_path)

    Z_pca, pca, pca_explained, pca_loadings = pca_reduce(X_scaled, feature_names, var_threshold)
    pca_load_path = os.path.join(figures_dir, "pca_loadings_heatmap.png")
    plot_heatmap(pca_loadings, feature_names, [f"PC{i+1}" for i in range(pca_loadings.shape[1])], "PCA主成分载荷热力图", pca_load_path)
    k_best_2, ks2, sc2 = silhouette_k_search(Z_pca, k_min, k_max, random_state, n_init)
    curve2 = os.path.join(figures_dir, "kmeans_on_pca_silhouette_curve.png")
    plot_silhouette_curve(ks2, sc2, "方式二：PCA降维后聚类的轮廓系数-簇数关系", curve2)
    km2 = KMeans(n_clusters=k_best_2, random_state=random_state, n_init=n_init)
    labels2 = km2.fit_predict(Z_pca)
    sil2_overall, sil2_cluster, sil2_samples = silhouette_details(Z_pca, labels2)
    scatter2 = os.path.join(figures_dir, "pca_scatter.png")
    plot_scatter_2d(Z_pca[:, :2], labels2, "PCA二维散点图（按簇着色）", scatter2, ("PC1", "PC2"))
    stability2 = stability_variance(Z_pca, k_best_2, [101, 123, 231, 345, 567], n_init)
    # 将PCA簇中心重构回原标准化特征空间以命名簇
    centers2_recon = km2.cluster_centers_ @ pca.components_
    clusters2_label = label_clusters_by_centers(centers2_recon, feature_names)
    # 命名主成分
    pca_named = name_all_dimensions(pca_loadings, feature_names, is_factor=False)
    t4 = time.perf_counter()

    # 4) 因子分析 + Varimax 后聚类
    t5 = time.perf_counter()
    Z_fa, fa_load_rot, fa_load_raw, n_factors, fa_explained = factor_analysis_reduce(X_scaled, feature_names, var_threshold)
    fa_load_path = os.path.join(figures_dir, "fa_loadings_heatmap.png")
    plot_heatmap(fa_load_rot, feature_names, [f"F{i+1}" for i in range(fa_load_rot.shape[1])], "因子载荷热力图（Varimax旋转）", fa_load_path)
    k_best_3, ks3, sc3 = silhouette_k_search(Z_fa, k_min, k_max, random_state, n_init)
    curve3 = os.path.join(figures_dir, "kmeans_on_fa_silhouette_curve.png")
    plot_silhouette_curve(ks3, sc3, "方式三：因子分析后聚类的轮廓系数-簇数关系", curve3)
    km3 = KMeans(n_clusters=k_best_3, random_state=random_state, n_init=n_init)
    labels3 = km3.fit_predict(Z_fa)
    sil3_overall, sil3_cluster, sil3_samples = silhouette_details(Z_fa, labels3)
    scatter3 = os.path.join(figures_dir, "fa_scatter.png")
    plot_scatter_2d(Z_fa[:, :2], labels3, "因子得分二维散点图（按簇着色）", scatter3, ("F1", "F2"))
    stability3 = stability_variance(Z_fa, k_best_3, [101, 123, 231, 345, 567], n_init)
    # 将FA簇中心重构回原标准化特征空间以命名簇（Z_fa @ loadings^T ≈ X_std）
    centers3_recon = km3.cluster_centers_ @ fa_load_rot.T
    clusters3_label = label_clusters_by_centers(centers3_recon, feature_names)
    # 命名因子
    fa_named = name_all_dimensions(fa_load_rot, feature_names, is_factor=True)
    t6 = time.perf_counter()

    # 5) 可解释性评分
    interp1 = interpretability_score_kmeans(centers1_z)
    interp2 = interpretability_score_from_loadings(pca_loadings)
    interp3 = interpretability_score_from_loadings(fa_load_rot)

    # 汇总与报告
    figures = {
        "kmeans_curve": curve1,
        "pca_scree": scree_path,
        "pca_loadings": pca_load_path,
        "pca_scatter": scatter2,
        "pca_curve": curve2,
        "fa_loadings": fa_load_path,
        "fa_scatter": scatter3,
        "fa_curve": curve3,
    }

    summary = {
        "kmeans": {
            "silhouette": sil1_overall,
            "runtime": (t2 - t1),
            "stability_var": stability1,
            "interpretability": float(interp1),
            "best_k": float(k_best_1),
        },
        "pca": {
            "silhouette": sil2_overall,
            "runtime": (t4 - t3),
            "stability_var": stability2,
            "interpretability": float(interp2),
        },
        "fa": {
            "silhouette": sil3_overall,
            "runtime": (t6 - t5),
            "stability_var": stability3,
            "interpretability": float(interp3),
        },
    }

    md_path = write_markdown_report(
        out_dir=out_dir,
        figures=figures,
        summary=summary,
        pca_info={
            "n_components": int(pca.n_components_),
            "explained_ratio_full": pca_full.explained_variance_ratio_,
        },
        fa_info={
            "n_factors": int(n_factors),
            "approx_explained": fa_explained,
        },
        code_params={
            "random_state": str(random_state),
            "n_init": str(n_init),
        },
        cluster_labels_direct=clusters1_label,
        cluster_labels_pca=clusters2_label,
        cluster_labels_fa=clusters3_label,
        pca_named=pca_named,
        fa_named=fa_named,
    )

    return md_path, figures


def build_argparser() -> argparse.ArgumentParser:
    """
    构建命令行参数解析器。
    """
    parser = argparse.ArgumentParser(
        description="美元期货指数聚类分析（三种方案：直接、PCA后、FA后）",
    )
    default_csv = os.path.join(os.path.dirname(__file__), "美元指数期货历史数据.csv")
    default_out = os.path.join(os.path.dirname(__file__), "outputs")
    parser.add_argument("--csv", type=str, default=default_csv, help="CSV数据路径")
    parser.add_argument("--out-dir", type=str, default=default_out, help="输出目录（含figures）")
    parser.add_argument("--k-min", type=int, default=2, help="K最小值")
    parser.add_argument("--k-max", type=int, default=10, help="K最大值")
    parser.add_argument("--random-state", type=int, default=123, help="随机种子（主流程固定为123）")
    parser.add_argument("--n-init", type=int, default=10, help="KMeans初始化次数")
    parser.add_argument("--var-threshold", type=float, default=0.7, help="PCA/FA累计解释度阈值")
    return parser


def main() -> None:
    """
    命令行入口：解析参数并运行完整流程，输出图表与报告。
    """
    set_chinese_font()
    parser = build_argparser()
    args = parser.parse_args()

    md_path, figures = run_pipeline(
        csv_path=args.csv,
        out_dir=args.out_dir,
        k_min=args.k_min,
        k_max=args.k_max,
        random_state=args.random_state,
        n_init=args.n_init,
        var_threshold=args.var_threshold,
    )

    print("报告生成：", md_path)
    for k, v in figures.items():
        print(f"图表[{k}]：{v}")


if __name__ == "__main__":
    main()